Skip to content

Configurable backend and per-solve scratch allocator - #467

Open
lkdvos wants to merge 25 commits into
mainfrom
reuse-buf
Open

Configurable backend and per-solve scratch allocator#467
lkdvos wants to merge 25 commits into
mainfrom
reuse-buf

Conversation

@lkdvos

@lkdvos lkdvos commented Jul 20, 2026

Copy link
Copy Markdown
Member

Every local update of an MPS algorithm contracts a handful of intermediate tensors that die before the next update.
This PR stops handing those to the garbage collector and serves them from a dedicated allocator instead, and while at it makes the contraction backend a setting on every algorithm.

Two settings, in two different places

backend is a field on the algorithm, defaulting to Defaults.backend().
It is stateless, and TensorOperations.DefaultBackend is a placeholder that resolves to an implementation from the types of the tensors involved, so it is already the right choice for a GPU-backed state.

The allocator is not a field. Which one is appropriate depends on two things the algorithm cannot know:

  • Where the tensors live. Host memory can be served by an allocator that bypasses Julia's memory manager; a device-backed state needs its scratch space on the device. This is only known once there is a state in hand.
  • Whether it will be shared between tasks. BufferAllocator is a bump buffer with a mutable offset, so it is only safe when a single task owns it. Concurrent work gets a ManualAllocator instead, which holds no state at all.

So MPSKit.default_allocator(state, scheduler) is asked once at the start of each solve, and the answer serves every local update of that solve.
This is what settles the concern from the first version of this PR: the algorithms stay thread-safe, and reusing one across a parallel parameter sweep is fine — there is no scratch space on the algorithm to share.
Anything MPSKit does not recognise falls back on DefaultAllocator, which allocates through the storage type itself and is therefore correct on any device.

Buffering can be turned off with MPSKit.Defaults.set_buffering!(false), which trades the speed back for lower memory use.
It is a compile-time preference, so Julia has to be restarted.

What it buys

Two-site DMRG ground state, spin-1 Heisenberg, L = 100, χ = 128, 5 sweeps, single Julia thread, MKL with 16 threads:

metric (steady state) DefaultAllocator BufferAllocator change
wall time 65.5 s 50.4 s −23%
GC time 24.2 s (37%) 10.3 s (20.5%) −57%
allocated 110.85 GiB 39.85 GiB −64%
GC pauses ~311 ~137 −56%
ground-state energy −138.940086126678 −138.940086126678 Δ = 0

Of the ~111 GiB a sweep churns through, ~71 GiB is recyclable effective-operator scratch.
The residual ~40 GiB is Krylov basis vectors, the gauge-step SVD, environment updates and result tensors, which this allocator does not manage and which is therefore the floor.

Smaller measurements from the same machine (@allocated, so read as ratios): finite TDVP 513 → 325 MiB, BUG 886 → 792 MiB, CBE-DMRG 1263 → 1170 MiB.
BUG gains least because much of its cost sits in the basis-augmentation step, which is a factorization rather than an effective-operator contraction.

Scope

Wired: DMRG/DMRG2, IDMRG/IDMRG2, TDVP/TDVP2, BUG, VUMPS, VOMPS, GradientGrassmann, DMRG3S, the bond-expansion algorithms, ChepigaAnsatz2, FiniteExcited, exact_diagonalization, dynamical DMRG, and the approximate/leading_boundary variants of the above.
calc_galerkin takes them too, so the Galerkin error no longer bypasses the buffer.

Not wired, and left for a follow-up: fidelity_susceptibility and the quasiparticle excitations, neither of which has an algorithm object carrying a backend today.

Notes for review

  • The effective-operator structs (MPODerivativeOperator, the JordanMPO_* operators, ProjectionDerivativeOperator) now carry the backend and allocator their matvec uses. Threading the 37 @plansor calls in the Jordan operators is what makes any of this pay: before that, the plumbing was inert for MPOHamiltonian and the two allocators were bit-identical in allocations.
  • The allocator is obtained at the entry point of a solve and threaded down. Obtaining it anywhere else silently loses the reuse — it went wrong three times during review (per IDMRG iteration, per site in changebond!, and not at all in calc_galerkin), each time without any test noticing, since the answers are unaffected.
  • Steps with no effective-operator contraction (a plain gauge!, RandExpand's changebond!) accept the allocator and drop it, so a sweep can pass it unconditionally.
  • Test coverage lives next to the algorithms it belongs to rather than in a file of its own: the finite DMRG/timestep testsets check that reusing one algorithm object reproduces the answer, and the infinite VUMPS, GradientGrassmann and TDVP testsets are parametrized over the serial and dynamic schedulers, since that is what picks the allocator. test/misc/allocator.jl covers the selector itself, and test/gpu/{cuda,amd}/allocator.jl covers only what a real device can settle — that a device-backed state is never handed a host allocator.

@lkdvos
lkdvos marked this pull request as draft July 20, 2026 19:09
@lkdvos
lkdvos requested review from borisdevos and leburgel July 20, 2026 19:09
@codecov

codecov Bot commented Jul 20, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 87.69231% with 40 lines in your changes missing coverage. Please review.

Files with missing lines Patch % Lines
src/algorithms/derivatives/mpo_derivatives.jl 64.44% 16 Missing ⚠️
.../algorithms/derivatives/hamiltonian_derivatives.jl 91.80% 5 Missing ⚠️
src/algorithms/grassmann.jl 66.66% 5 Missing ⚠️
...c/algorithms/derivatives/projection_derivatives.jl 73.33% 4 Missing ⚠️
src/utility/defaults.jl 20.00% 4 Missing ⚠️
src/algorithms/approximate/vomps.jl 72.72% 3 Missing ⚠️
src/algorithms/excitation/chepigaansatz.jl 80.00% 2 Missing ⚠️
src/algorithms/changebonds/randexpand.jl 50.00% 1 Missing ⚠️
Files with missing lines Coverage Δ
src/MPSKit.jl 100.00% <ø> (ø)
src/algorithms/ED.jl 100.00% <100.00%> (ø)
src/algorithms/approximate/fvomps.jl 95.83% <100.00%> (+4.52%) ⬆️
src/algorithms/approximate/idmrg.jl 96.77% <100.00%> (+0.05%) ⬆️
src/algorithms/changebonds/optimalexpand.jl 98.87% <100.00%> (+0.02%) ⬆️
src/algorithms/changebonds/sketchedexpand.jl 98.43% <100.00%> (ø)
src/algorithms/changebonds/vumpssvd.jl 100.00% <100.00%> (ø)
src/algorithms/derivatives/derivatives.jl 92.15% <100.00%> (+0.49%) ⬆️
src/algorithms/excitation/dmrgexcitation.jl 100.00% <100.00%> (ø)
src/algorithms/groundstate/dmrg.jl 92.85% <100.00%> (+0.22%) ⬆️
... and 22 more

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@lkdvos
lkdvos marked this pull request as ready for review July 20, 2026 22:42
@lkdvos

lkdvos commented Jul 21, 2026

Copy link
Copy Markdown
Member Author

Benchmark: DefaultAllocator vs BufferAllocator

End-to-end DMRG2 ground-state run measuring the effect of reusing a BufferAllocator across the sweep, driven by benchmark/allocator_gc.sbatch / benchmark/allocator_gc_run.jl.

Setup: spin-1 Heisenberg chain (Trivial), L=100, χ=128, 5 fixed sweeps, krylovdim=20, non-adaptive eigensolver (identical work for both allocators). Single Julia thread (so GC is single-threaded and its wall time is a clean number), MKL with 16 threads. julia 1.12.6.

metric (steady-state) DefaultAllocator BufferAllocator change
wall time 65.5 s 50.4 s −23%
GC time 24.2 s (37%) 10.3 s (20.5%) −57%
allocated 110.85 GiB 39.85 GiB −64% (−71 GiB)
GC pauses ~311 ~137 −56%
full collections 22 9 −59%
ground-state energy −138.940086126678 −138.940086126678 Δ = 0

Of the ~111 GiB a sweep churns through, ~71 GiB is recyclable effective-operator contraction scratch — with BufferAllocator it now comes from the reused buffer instead of the GC heap. The residual ~40 GiB (Krylov basis vectors, the gauge-step SVD, environment updates, result tensors) is not managed by the operator allocator and is therefore identical for both — that's the floor.

Raw output
model=Heisenberg S=1  L=100  χ=128  sweeps=5  krylovdim=20  reps=3
julia=1.12.6  julia_threads=1  gc_threads=1  BLAS_threads=16
DefaultAllocator | rep 1/3 |    66.79 s | GC   25.48 s ( 38.1%) |   110.85 GiB |    320 pauses |  23 full
DefaultAllocator | rep 2/3 |    65.62 s | GC   24.24 s ( 36.9%) |   110.85 GiB |    310 pauses |  22 full
DefaultAllocator | rep 3/3 |    65.41 s | GC   24.14 s ( 36.9%) |   110.85 GiB |    313 pauses |  22 full
BufferAllocator  | rep 1/3 |    50.91 s | GC   10.34 s ( 20.3%) |    39.85 GiB |    135 pauses |   9 full
BufferAllocator  | rep 2/3 |    50.37 s | GC   10.33 s ( 20.5%) |    39.85 GiB |    137 pauses |   9 full
BufferAllocator  | rep 3/3 |    50.38 s | GC   10.32 s ( 20.5%) |    39.85 GiB |    137 pauses |   9 full
energy check: |Δ| = 0.00e+00

Comment thread src/algorithms/derivatives/hamiltonian_derivatives.jl
Comment thread src/algorithms/derivatives/hamiltonian_derivatives.jl Outdated

@leburgel leburgel left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks great to me. Aside from the comments I am quite curious: do you have any idea of how this affects the compilation time?

Comment thread src/algorithms/groundstate/dmrg.jl Outdated
Comment thread src/algorithms/derivatives/hamiltonian_derivatives.jl

@lkdvos lkdvos left a comment

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure this affects compilation times that much, in the sense that of course we have to compile stuff multiple times for different backends and allocators, but since that was previously not possible this isn't really a regression.

I'm still playing around a bit with some compilation time stuff too though, for example I think the higher-level algorithm implementations don't actually need to be specialized which might recover some of this.

Comment thread src/algorithms/derivatives/hamiltonian_derivatives.jl
@lkdvos

lkdvos commented Jul 30, 2026

Copy link
Copy Markdown
Member Author

Update here: I'm still working on what I want to do with the design here, and for now I'm thinking that this might be overengineered right now, so I'm going to try and slim this down slightly.

@lkdvos

lkdvos commented Jul 31, 2026

Copy link
Copy Markdown
Member Author

Alright, update again here, I think this is now more or less settling on a design I'm happy with. I kept the backends selectable, being a hook for making changes, but actually decided against making the allocator a configurable thing, mostly because it is just too easy to accidentally do something wrong with that. I think it is fair to always take the bumper where possible, and this now also handles the GPU case, so the main thing left for me to do is clean up the tests and rebase. Also note that this now means that parameter sweeps in parallel would work perfectly fine, and would probably have way better performance because the GC needs to run less often.
I'll try to do some more before/after tests to actually see if this matters.

Comment thread src/algorithms/groundstate/idmrg.jl Outdated
Comment thread src/algorithms/groundstate/dmrg.jl Outdated
Comment thread src/algorithms/derivatives/projection_derivatives.jl
Comment thread src/algorithms/ED.jl Outdated
lkdvos and others added 14 commits August 3, 2026 18:18
`Defaults.backend()` names the default contraction backend in one place, so the
algorithms that grow a `backend` field can default it consistently.
`TensorOperations.DefaultBackend` is a placeholder rather than a CPU backend --
the implementation is selected from the types of the tensors involved -- so this
is already the right choice for GPU-backed states.

`Defaults.buffering` / `set_buffering!` is the opt-out for the dedicated scratch
space added next, for when memory rather than time is the binding constraint.
It is a compile-time preference, stored in `LocalPreferences.toml`.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The local updates of an MPS algorithm allocate a fair number of intermediates,
all of which die before the next update. Serving them from an allocator rather
than from Julia's memory manager keeps them out of the garbage collector's way.

Which allocator is appropriate depends on where the tensors live and on whether
more than one task will be sharing it, neither of which is a property of the
algorithm: `default_allocator(x, scheduler)` is therefore asked at the start of a
solve, and the answer used for all of its local updates. Host memory gets a
`BufferAllocator` when a single task owns it and a `ManualAllocator` when it is
shared; any other storage type falls back on `DefaultAllocator`, which allocates
through the storage type itself and is therefore correct on any device.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The effective local operators now carry the backend and allocator their matvec
should use, and their constructors and `prepare_operator!!` accept them. This is
what lets a single reused buffer capture all of the intermediates of the
effective-Hamiltonian applications across an entire sweep.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The bond-expansion algorithms grow a `backend` field and thread it, together with
an allocator, through their contractions. `changebond!` and `gauge!` take them as
keyword arguments, so a `DMRG`/`TDVP` sweep can hand its own down into the
expansion and gauge steps rather than each step picking its own.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`DMRG`/`DMRG2` grow a `backend` field, and the sweep asks `default_allocator` once
for a serial allocator that then serves every local update: a single buffer
captures the effective-operator intermediates of the whole sweep. On a spin-1
Heisenberg chain (L=100, chi=128) this cuts allocations by 64%, garbage-collection
time by 57% and wall time by 23%, at an unchanged ground-state energy.

`IDMRG`/`IDMRG2` and the `approximate` variants get the same treatment.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Both integrators grow a `backend` field. The finite sweeps take a single serial
allocator for all of their local updates; for the infinite `TDVP` the scheduler is
read once at the entry point, so that the allocator it selects -- a buffer when
serial, a manual allocator when the sites are spawned concurrently -- is
inferable, and the AC and C sweeps share it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Both grow a `backend` field, and their unit-cell sweeps take the allocator
matching the scheduler they spawn with.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`fg` takes the backend and the scheduler as keyword arguments rather than reading
the scheduler itself, so the allocator it selects matches the concurrency of the
gradient it is computing.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`ChepigaAnsatz2`, `FiniteExcited`, `exact_diagonalization` and the dynamical-DMRG
propagator grow a configurable backend and draw their local-update scratch space
from an allocator, as the other algorithms now do.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Document why the allocator is asked for per solve rather than configured on the
algorithm: it depends on where the tensors live and on whether it will be shared
between tasks, and neither is known to the algorithm. Also note that `backend`
needs no such treatment, since `DefaultBackend` already resolves to cuTENSOR for
`CuArray`-backed states.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…orithms

`test/misc/allocator.jl` covers `default_allocator` itself: the host/device split,
the serial-vs-shared choice and the buffering preference. `test/algorithms/allocator.jl`
checks that the algorithms give the same answers with and without a buffer, and
`test/gpu/*/algorithms.jl` runs a shared set of algorithms on device-backed states,
which is what pins down that the device path stays on `DefaultAllocator`.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Same treatment as the other integrators: a `backend` field, and a single serial
allocator obtained once per `timestep!` that serves the effective-Hamiltonian
application of every local update. The bond cut and the basis augmentation are
factorizations rather than effective-operator contractions, so -- as with a plain
`gauge!` -- there is nothing to route to them.

Measured on a spin-1 Heisenberg chain (L=20, chi=48), one BUG step allocates
792 MiB with the buffer against 886 MiB without it; the smaller relative gain than
TDVP's (325 against 513 MiB) is because BUG spends more of its time in the
augmentation step, which the allocator does not serve.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ctions

`ProjectionDerivativeOperator` was the one derivative operator whose matvec still
contracted with a bare `@plansor`, so the Galerkin error and the `FiniteExcited`
eigensolves left their intermediates to the garbage collector even when the sweep
was holding a buffer. It now carries the backend and allocator like the others do,
and the constructors pass them in.

Reported by @borisdevos in review of #467.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
lkdvos and others added 6 commits August 3, 2026 18:49
The Galerkin error is one full effective-operator application per site, which was
going through the default backend and the garbage collector. It now takes both as
keyword arguments and forwards them to `AC_projection`; DMRG passes the allocator
its sweep already holds, and the infinite algorithms pass their backend (their
convergence check runs outside the sweep that owns the allocator).

Reported by @borisdevos in review of #467.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`localupdate_step!` asked `default_allocator` on every outer iteration, so each
sweep grew a fresh buffer from empty and dropped it again -- the opposite of the
reuse the buffer exists for. The allocator now lives on `IDMRGState`, next to the
timer output, obtained once when the solve starts and carried along; the initial
Galerkin error uses it too.

Reported by @borisdevos in review of #467.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`Bk <: AbstractBackend`, as `MPODerivativeOperator` and `PrecomputedDerivative`
already have it. This buys nothing from the compiler -- Julia specializes on the
field type either way -- but it does give an explicit error rather than a late
`MethodError` if something that is not a backend ends up in that field.

Suggested by @borisdevos in review of #467.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`test/setup/gpu_algorithms.jl` reached for `TensorOperations` directly, which is not a
test dependency, so the CUDA/AMDGPU algorithm tests errored on load. Take the
allocator types from MPSKit, as `test/misc/allocator.jl` already does.

While here, run the device timestep for `BUG` as well as `TDVP`. Verified on an
RTX A6000: 19/19 pass.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`changebond!` obtained a fresh allocator on every call, i.e. once per site: with
CBE the buffer was rebuilt and thrown away at every local update, which is the
opposite of what it is for. It now takes the allocator as a keyword argument,
defaulting to one of its own so that a standalone `changebonds!` still works, and
the `DMRG`/`TDVP` sweeps pass the one they are already holding. `RandExpand` has no
effective-operator contraction to serve, so it accepts and drops it, as a plain
`gauge!` does.

Measured on a spin-1 Heisenberg chain (L=30, expand 8, truncate to 48), two
CBE-DMRG sweeps allocate 1170 MiB against 1263 MiB before.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
lkdvos and others added 4 commits August 4, 2026 09:05
… the operator

Every operator with a `prepare_operator!!` method now stores its own backend and
allocator, so passing them in again was redundant -- and every call site passed
exactly the pair the operator had just been constructed with. Both
`prepare_operator!!` and `prepared_operator_type` lose the extra arguments; the
prepared type's `B`/`A` parameters come from the unprepared operator's.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`test/algorithms/allocator.jl` re-ran DMRG, DMRG2, TDVP, BUG, CBE, IDMRG and
changebonds that the algorithm suites already cover. Only two claims were actually
new, and they now live next to the tests they belong to:

- re-using one algorithm object across solves reproduces the answer (it carries no
  scratch space of its own), in the finite `DMRG` and `timestep` testsets;
- the scheduler choice -- which is what decides between a buffer and a shared
  manual allocator -- changes no number, by parametrizing the existing infinite
  `VUMPS`, `GradientGrassmann` and `Infinite TDVP` testsets over both.

`with_scheduler` and the serial/dynamic pair move to `TestSetup`, so there is one
failure-safe implementation instead of the ad-hoc `set_scheduler!` save/restore
`approximate.jl` was doing. `MPSKit.Defaults.scheduler` is a `Ref` rather than a
`ScopedValue`, hence the explicit restore.

The hand-rolled inference helpers are gone in favour of `@testinferred`; the one
remaining `return_type` call is the abstract-signature case, which no call-based
macro can express, and matches how `src/` already spells it.

The GPU file keeps only what a device can prove -- that a device-backed state is
never handed a host allocator -- and drops the `find_groundstate`/`timestep` runs,
which duplicated CPU coverage. Renamed to `allocator.jl` accordingly, since it no
longer tests algorithms.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…testsets

With two loop variables the `@testset for` body is a single scope for all
iterations, so binding the starting state to `ψ` fed the previous iteration's
result back into `repeat(ψ, unit_cell_size)`: the second scheduler arm at unit
cell 3 got a length-9 state against a length-3 Hamiltonian ("lengths must match").

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@lkdvos
lkdvos marked this pull request as ready for review August 4, 2026 17:08
@lkdvos
lkdvos requested review from borisdevos and leburgel August 4, 2026 17:09
@lkdvos

lkdvos commented Aug 4, 2026

Copy link
Copy Markdown
Member Author

I think this is now actually in a merge-able state. I'm kind of happy with taking the allocators out of the hands of configuration, and it seems like it doesn't complicate the code too much, so I think this is acceptable

@lkdvos lkdvos changed the title use backends and allocators in DMRG Configurable backend and per-solve scratch allocator Aug 4, 2026
The testsets in `InfiniteMPS ground state` share one mutable `ψ`: each rebinds it,
and its `unit_cell_size == 1` iteration resets it to a fresh length-1 state, which
is what makes the `3` iteration's `repeat(ψ, 3)` correct. Parametrizing two of them
over the scheduler broke that convention -- they no longer rebind `ψ`, so they read
whatever a neighbouring testset had left there, and `repeat`ing an already-repeated
state gave a length-9 state against a length-3 Hamiltonian.

They now build their own starting state and touch the shared `ψ` not at all.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants